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dient method (modified implicit midpoint rule) where we succeeded 
to preserve essential geometric properties and the final results have 
a relatively simple form. In the case of one-dimensional Hamiltonian 
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1 Introduction 



The motivation for introducing "locally exact" discretizations is quite natu- 
ral. Considering a numerical scheme for a dynamical system with periodic 
solutions (for instance: a nonlinear pendulum) we may ask whether the nu- 
merical scheme recovers the period of small oscillations. For a fixed finite 
time step the answer is usualy negative. However, many numerical scheme 
(perhaps all of them) admit modifications which preserve the period of small 
oscillations for a fixed (not necessarily small) time step h. An unusual feature 
of our approach is that instead of taking the limit — > 0, we consider the 
limit X (where x is the stable equilibrium). The next step is to consider 
a linearization around any fixed x (then, in order to preserve the condition 
Xn ~ some evolution of x is necessary). In other words, we combine 
two known procedures: the approximation of nonlinear systems by linear 
equations and explicit exact discretizations of linear equations with constant 
coefficients. An important novelty consists in applying these procedures to a 
modified numerical scheme containing free functional parameters. The case 
of small oscillations was presented in [8]. More results for one-dimensional 
Hamiltonian systems can be found in [9l [11] (by one-dimensional Hamilto- 
nian system we mean a Hamiltonian system with one degree of freedom). 
The results are very promising. It seems that a new, very accurate method 
is emerging. In this paper we extend our approach on the case of multi- 
dimensional canonical Hamiltonian systems. Our method works perfectly 
for discrete gradient schemes. We succeeded to modify the discrete gradient 
scheme in a locally exact way without spoiling its main geometric property: 
the exact conservation of the energy integral. It is well known that preser- 
vation of geometric properties by numerical algorithms is of considerable 
advantage [151 [H] • 

A notion identical with our local exactness has been proposed a long time 
ago [29], see also [22]. Recently, similar concept appeared under the name 
of linearization-preserving preprocessing [25J, see also below (section 13. 2p . 
Our approach has also some similarities with the Mickens approach [2B] . 
Gautschi-type methods [121 [IZ] and, most of all, with the exponential in- 
tegrators technique [3l [161 [2I|- The definition of exponential integrators is 
so wide (e.g., "a numerical method which involves an exponential function 
of the Jacobian", [16]) that our numerical schemes can be considered as 
special exponential integrators. In spite of some similarities and overlap- 
pings, our approach differs from all methods mentioned above. In particu- 
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lar, according to our best knowledge, discrete gradient schemes have never 
been treated or modified in the framework of exponential integrators and/or 
linearization-preserving preprocessing. Main new results of this paper consist 
in constructing energy-preserving locally exact discrete gradient schemes for 
arbitrary multidimensional Hamiltonian systems, see section |6l 



2 Exact discretization of linear systems 

We consider an ordinary differential equation (ODE) with the general so- 
lution x(t) (satisfying the initial condition xito) = Xq), and a difference 
equation with the general solution x„. The difference equation is the exact 
discretization of the considered ODE if Xn = x{tn). 

It is well known that any linear ODE with constant coefficients admits the 
exact discretization in an explicit form [30] , see also [H [TJ [26] . We summarize 
these results as follows, compare [11] (Theorem 3.1). 

Proposition 2.1. Any linear equation with constant coefficients, represented 
in the matrix form by 

dx 

— = Ax + b, (2.1) 
dt 

(where (t) e W^, b = const e and A is a constant invertible real 

d X d matrix) admits the exact discretization given by 

Xn+i -Xn = (e'^"^ - 1)A-^ {AXn + 6) , (2.2) 

where hn = tn+i — is the time step and 1 is the identity matrix. 

Corollary 2.2. We may look at the exact discretization ^2. 31) as a modifica- 
tion of the forward Euler scheme. Indeed, we can rewrite ^2.2) as 

A;;^ {Xn+l - Xn) = AXn + 6 , (2.3) 

where An, defined by 

A„ = A-i(e'^"^-l) , (2.4) 
is a matrix "perturbation" of the time step hn. Indeed, An = hnl + 0{hD. 
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As a special case we consider the multidimensional harmonic oscillator 
equation driven by a constant force: 

+ fi^x = a , (2.5) 

where x = x{t) G M*" and f2 is a given invertible constant m x m matrix 
and a = const G M™. It is convenient to represent f l2.5p as the following first 
order system 

x=p, p = —il'^x + a . (2.6) 

Proposition 2.3 ([6], Prop. 19). The exact discretization of the system \2.0\) 
is given by 



COS flh„. 0, ^ sin Qhn \ f Xr, \ . f 2Q sin^ ^ a 

(2.7) 



Pn+i ) V sin cosVthn J \ Pn J ^ \ ^ ^sinf2/i„a 



(2.8) 



where hn is an arbitrary variable time step. 

Corollary 2.4 (|6j, Prop. 20). The formulas l[2. 7| ) can be rewritten in the 
following equivalent form: 

K\Xn+l -Xn) = \ [Pn+l + Pn) , 

- Pn) = -l^'^iXn+1 + Xn) + tt . 

where 

5n = 2i7~^ tan ^ = hntanc ^ (2-9) 

Here and below we use notation tanc(z) = tanz. It is worthwhile to 
notice that tanc(2;) is an even function and depends analytically on z"^. 

As one could expect, the exact discrete harmonic oscillator equations 
preserve exactly the total energy. 

Proposition 2.5 ([6J, Prop. 22). // Q is a symmetric matrix (Q^ = 
then 

1 1 

In ■=2^n\ Pn) + l^M ^'^Xn) - (x„ | o) (2.10) 

is an integral of motion (i.e., In+i = In) of the discrete multidimensional har- 
monic oscillator equations Ii2.8\) . Here the bracket denotes the scalar product 
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Exact discretization seems to be of limited value because, in order to apply 
it, we need to know the explicit solution of the considered system. However, 
there exist non-trivial applications of exact discretizations. In the case of the 
classical Kepler problem we succeeded in using the exact discretization of the 
harmonic oscillator in two different ways, obtaining two different numerical 
integrators preserving all trajectories and integrals of motion [H O E]- The 
exact discretization of the harmonic oscillator equation can also be used to 
in the integration of some partial differential equations by Fourier transfor- 
mation [6j. 



3 Locally exact numerical schemes 

The central topic of our paper is another fruitful direction of using exact 
integrators, namely the so called locally exact discretizations [6l [9] . 



3.1 Motivating example 

First, we recall our earlier results concerning one-dimensional Hamiltonian 
systems p = —V'{x), x = p, see [H E]- We tested the following class of 
numerical integrators 

7 = o {Pn+1 + Pn) ■ 

(3.1) 

Pn+1 - Pn _ V{Xn+l) - VjXn) 

where 6n is a function defined by 

2 hnUln 



5n = — tan^ , a;„ = ^V"{x) , (3.2) 

and, in general, x may depend on n. For simplicity, we formally assume 
V"{x) > 0. However, in the case of non-positive V"{x) one can use the same 
formula, for details and final results see [S]. The simplest choice is x = Xq, 
where V'{xq) = (small oscillations around the stable equilibrium). In this 
case Sn does not depend on n. The resulting scheme, known as MOD-GR, 
was first presented in [8] . In [9] we considered the (GR-LEX) and 

its symmetric (time- reversible) modification x = |(a;„ + x^+i) (GR-SLEX). 
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In both cases x is changed at every step. The names MOD-GR, GR-LEX, 
GR-SLEX first appeared in [IHl EH. 

Any numerical scheme from the family defined by (13.11) and (13. 2p is lo- 
cally exact (or, more precisely, locally exact at x), which means that its 
linearization at x coincides with the exact discretization of the linearized 
system. Indeed, the linearized system, given hj p = —V'{x) — V"{x)C,, 
C, = p, is a. particular case of (12. 6p and admits the exact discretization (I2.8p . 
(12. 9p . where we have to identify Q'^ = V"{x), a = —V'{x). This exact dis- 
cretization coincides with the linearization of (13. ip . obtained by substituting 
V{xn) ~ V{x) + V'{x)^n + ^V'{^)^ny Compare [9]. In this paper we are going 
to present a generalization of this approach. 

3.2 Local exactness 

Motivated by the results of [H |9] we propose the following definitions. 

Definition 3.1. A numerical scheme Xn+i = ^{Xn,hn) for an autonomous 
equation X = F{x) is locally exact atx if its linearization around x is identical 
with the exact discretization of the differential equation linearized around x. 

The scheme MOD-GR is locally exact at a stable equilibrium only, GR- 
LEX is locally exact at Xn (for any n), and, finally, GR-SLEX is locally exact 
at ^{xn + Xn+i) (for any n). In the case of implicit numerical schemes the 
function $(a:„, /i„) is, of course, an implicit function. 

Definition 3.2. A numerical scheme Xn+i = ^{Xn,hn) for an autonomous 
equation x = F{x) is locally exact if there exist a sequence Xn such that 
Xn — Xn = 0{hn) and the scheme is locally exact at Xn (for any n). 

Therefore GR-LEX and GR-SLEX are locally exact. Similar concept 
("linearization-preserving" schemes) has been recently formulated in [25] 
(Definition 2.1). An integrator is said to be linearization-preserving if it 
is linearization preserving at all fixed points. This definition is weaker than 
our Definition 13.21 All locally exact schemes are linearization-preserving 
but, in general, linearization-preserving scheme has not to be locally exact in 
our sense. For instance, the scheme MOD-GR (see [HI E]) is linearization- 
preserving (provided that V{x) has only one stable equilibrium) but is not 
locally exact. 
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3.3 Exact discretization of linearized equations 

As an immediate corollary from Definition 13.21 we have that the exact dis- 
cretization of the linearization of a given nonlinear system is locally exact. 
We confine ourselves to autonomous systems of the form 

X = F{x) (3.3) 
where x{t) G M^. If X is near a fixed x, then (13. 3p can be approximated by 

e = F'(x)e + F{x) (3.4) 
where F' is the Frechet derivative (Jacobi matrix) of F and 

$ = x-x. (3.5) 
The exact discretization of the approximated equation (13.41) is given by 

Ui = e'-'^'^'^L + (e'"^'(") - l) {F'{x)r' F{x) , (3.6) 

provided that det F'{x) ^ 0. We make this assumption here and throughout 
the paper. 

Corollary 3.3. The exact discretization of the linearization of x = F{x) 
around x is given by: 

Xn+i -x = e''"^'^^) {xn -x)+ (^e'^"^'(^) - l) {F'{x))~^ F{x) . (3.7) 



The scheme ^3. 7|) is locally exact 



If we put X = Xn into (13.71) (thus changing x at each step), then ^„ = 
and we obtain: 

Xn+l =Xn+ (e'^"^'(-") - l) {F'{Xn))-'F{Xn) (3.8) 

This method is well known as the exponential difference equation (see [29] , 
formula (3.4)) or, more recently, as the exponential Euler method 



+ hn Lpi{hnF' {Xn)) F{Xn) , (3.9) 



where ipi{z) = z ^{e^ — 1) (and V5i(0) := 1). Another possibihty, x = x^, 
where 

X+ = ^ (Xn+l+Xn) , (3.10) 

leads to 

-x„ = 2(F'(x+))-4anh Q/i„F'(^+)^ F(S+) . (3.11) 
A third natural possibility, x = Xn+i, yields the scheme: 

Xn+l = Xn + (l - e-'^-^'C-^+i)) {F'{Xn+l))''F{Xn+l) , (3.12) 

which may be rewritten in terms of (fi as 

Xn+l =Xn + K ipi{-KF' iXn+l)) F{Xn+l) , (3.13) 

Example 3.4. In the case p = —V'{x), x = p, we have 

^(^) = ( -n.) ) • ^'w = ( -A.) I ) • P'"' 

and the scheme ^3. 8\) assumes the form 

_ sm{hnUJn) _ 1 - cos(/ina;n) , 

sin(/i„w„) 

Pn+i = p„cos(/i„cu„) K (x„) , 



(3.15) 



where Un = \/V"{xn)- Global properties of this scheme are rather poor but 
it can serve as a very good predictor 

3.4 Linear stability of locally exact integrators 

Locally exact numerical schemes have excellent qualitative behavior around 
all fixed points of the considered system. 

Proposition 3.5. If the equation x = F{x) has a fixed point at x = x, then 
all its locally exact discretizations have a fixed point at x„ = x, as well. 
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Proof: If F{x) = 0, then equation (|3.7p becomes 

x„+i-x = e'^"^'(^)(x„-x) . (3.16) 

We require that the scheme x„+i = <I>(x„,/i) is a locahy exact discretization of 
F{x) = 0, i.e., its hnearization, given by 

Xn+l = ^{x, h) + $'(x, h){xn - x) , (3.17) 

coincides with ()3.16p . Hence 

^h^F'{x) ^ ^/(^^ ^ ^(^^ h)=x , (3.18) 

which means that S is a fixed point of the system Xn+i = ^{xn, h). □ 

Stability of numerical integrators can be roughly defined as follows: "the 
numerical solution provided by a stable numerical integrator does not tend to 
infinity when the exact solution is bounded" (see |3], p. 358). The integrator 
which is stable when applied to linear equations is said to be linearly stable. 
We may use the notion of A-stabihty (see, e.g., [H]): an integrator is said 
to be A-stable, if discretizations of stable linear equations are stable as well 
(equation i; = Ax is stable if ReA < 0). 

Making one more assumption (quite natural, in fact) that the discretiza- 
tion of a linear system is linear, we see that locally exact integrators are 
linearly stable. Indeed, solutions of any locally exact discretization have the 
same trajectories as corresponding exact solutions. 

Corollary 3.6. Any locally exact numerical scheme is linearly stable and, in 
particular, A-stable. 

Some simple examples illustrating this corollary will be given in sec- 
tion 14.51 In fact, much stronger result holds: a locally exact discretization 
yields the best (exact) simulation of a linear equation in the neighborhood of 
a fixed point. In particular, locally exact discretizations preserve any qual- 
itative features of trajectories of linear equations (up to round-off errors, of 
course) . 

Numerical experiments show that locally exact schemes are exceptionally 
stable also for some simple nonlinear systems [H [U [11], but we have no 
theoretical results concerning the stability (e.g., algebraic stability [19]) in 
the nonlinear case. 
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4 Locally exact modifications of popular one- 
step numerical schemes 

We are going to use local exactness as a criterion to select numerical schemes 
of higher accuracy from a family of non-standard integrators. Our working 
algorithm to derive such "locally exact modifications" of numerical integra- 
tors of the form (15. 2p assumes that (5„ depends only on x, p (or, in more 
general case, on x„) and hn. We usually consider only three different se- 
quences Xn, namely Xn = Xn, Xn = Xn+i and Xn = x'^ = |(x„ The 
problem of finding the best sequence for a given numerical scheme seems to 
be interesting but has not been considered yet. 

In this section we apply this procedure to a number of standard numerical 
methods. We do not claim that the resulting numerical schemes (sometimes 
new, sometimes already known) are much better than the corresponding 
schemes before the modification. Their accuracy should be considerably bet- 
ter but the computing cost is surely much higher. For every scheme separately 
one has to check the final outcome and to decide whether the modification 
is really profitable. Numerical simulations of small oscillations around stable 
fixed points and the results of [H [9] are quite promising. 

4.1 Locally exact explicit Euler scheme 

The explicit Euler scheme for (13. 3 p is given by Xn+i = x„ + hnF{Xn)- We 
postulate the following generalization of this scheme: 

Xn+l =Xn+ 5nF{Xn) (4.1) 

where 5„ is a variable matrix which is defined by the requirement that the 
linearization of (14. ip coincides with (13. 6p . We linearize (14. ip by substituting 
(ESD (i.e., Xn+x + i^): 

Ui=^n + ^n{F{x)+F\x)U (4.2) 

and assume that 5„ depends only on x and hn- Equation 14.21 is identical with 
f l33|) iff 

1 + 5nF'{x) = e'^"^'^^), 5n = (^e'*"^'(^) - l) {F'{x))-^ . (4.3) 

Once can easily see that both equations (14. 3 p are equivalent provided that 
F'{x) is non-degenerate. 
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Substituting 6n from fl4.3l) into fl4.ip we obtain a class of locally exact 
modifications of the explicit Euler scheme 

= x„ + (e^"^'(^) - 1) {F'{x)r'F{x„) . (4.4) 

It seems that the choice x = Xn is most natural in that case (then, by the 
way, fl4.4p coincides with (13. Sp ). 

Corollary 4.1. Locally exact modification of the explicit Euler scheme is 
given by 

= Xn + (e'^"^'(-") - 1) {F'iXn))-'F{Xn) (4.5) 

4.2 Locally exact implicit Euler scheme 

The implicit Euler scheme for (13. 3 p is given by Xn+i = Xn + hnF{Xn+i)- We 
postulate the following generalization of this scheme: 

Xn+l =Xn+ 5nF{Xn+l) (4.6) 

where 5n is a variable matrix which is defined by the requirement that the 
linearization of (14. 6 p coincides with (13. 6p . We linearize (14. 6 p by substituting 
dSSD (i.e.,a:„+x + e„): 

ia+l =ia+Sn {F{x) + F' {x)U,) , (4.7) 

i.e., (l — 5nF'{x))^n+i =in + ^nF{x), which is identical with (13. 6p iff 
(1-5„F'(S))"' = e^»^'(^) 

(4.8) 

{l~5^F'{x))-'5^ = (e''"^'^^) - l)iF'ix))-\ 

We easily see that 5„ can be computed independently from any of these two 
equations. Fortunately, the result is the same: 

Sn={l- e-^"^'(^)) {F'{x))-^ (4.9) 
Substituting Sn from (14.81) into (14. 6 p we get 

x„+i=a:„+(l-e-^"^'(^))(F'(^))-iF(a:„+i) , (4.10) 
or, in an equivalent way, 

Xn+l =Xn + hn (fi{-hnF' (x)) F{Xn+l) ■ (4.11) 

In this case it is not clear what is the most natural identification. We can 
choose either x = a;„, or x = Xn+i- 
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Corollary 4.2. Locally exact modification of the implicit Euler scheme is 
given either by 

=a;„+ (1 -e-'^"^'(^"))(F'(a:„))-iF(x„+i) , (4.12) 

or, by 

Xn+i = + (1 - e-'^-^'^^-^^)) {F'{xn+,))-'F{xn+i) • (4.13) 

4.3 Locally exact implicit midpoint rule 

We postulate the following extension of the implicit midpoint rule: 

x^^,-x^=6^F 1' ^-+^ +^" ') , (4.14) 



where Sn depends, as usual, on x and hn- The scheme fl4.14p is locally exact 
for 

6n = 2iF'ix))-Hanh( ^''^yA , (4.15) 



which will be shown, in more general framework, in section 14.61 In some 
applications the Taylor series may be useful: 

= 1 - ^hl{F'{x)r + ^ht,{F'{x)r + . . . (4.16) 
Thus locally exact modification of the implicit midpoint rule reads 

Xn+1 -Xn = hn{ tauhc 1 F ( 1 (4.17) 

The most natural choice of x seems to be at the midpoint, x = \ (Xn + Xn+i)- 
However, in order to diminish the computation cost, the choice X — Xyi can 
also be considered (because then the Jacobian F'{x) is evaluated outside 
iteration loops). 
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4.4 Locally exact trapezoidal rule 

We postulate the following extension of the trapezoidal rule: 

Xn+1-Xn=dn , (4.18) 

Linearizing fl4.18p around x we get 



(4.19) 



In section 14.61 we will show that the scheme fl4.18p is locally exact for 

Sn = 2{F'{x))-Heinh(^^^^^^] . (4.20) 



2 

Thus locally exact modification of the trapezoidal rule is given by 

/ hr.F'{x) \ F{Xn+l) + F{Xn) , . 

Xn+i -Xn = K\ tanhc j . (4.21) 

There are two natural choices of x. Either (in order to minimize the compu- 
tational costs) we can take x = Xn-, oi (in order to obtain a time- reversible 
scheme) we can take x = ^ (Xn + Xn+i)- 

4.5 A-stability of locally exact modifications 



In order to illustrate general results of section [331 we will apply four schemes 
presented above to the one-dimensional linear equation x = Xx. The locally 
exact explicit Euler scheme yields 

-Xn= (e'^"^ - 1) x„ . (4.22) 

The locally exact implicit Euler scheme yields 

Xn+l -Xn={l- e-'^"^ - 1) Xn+1 , (4.23) 

Both resulting equations are identical: Xn+i = e^^^Xn and yield the exact 
discretization. The implicit midpoint and trapezoidal rules yield an identical 
equation, namely 



l^tanh {xn+i + Xn) . (4.24) 

Computing Xn+i from (14.241) we get the exact discretization = e'^^^Xn 
again. In particular, if ReA < 0, then x„ — for n — )■ oo (for any constant 
time step hn = h = const). A-stability is evident in all these cases. 
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4.6 A large class of locally exact integrators 

All numerical schemes presented above are particular cases of the following 
class of numerical schemes for the equation x = F{x): 

Xn+l-Xn=S{x)'^{Xn,Xn+l) . (4.25) 

We assume the consistency conditions: 

^(rr,x) = F{x) , F'{x) = '^i{x,x) + '^2{x,x) , (4.26) 

where \l/i,\l/2 are partial Frechet derivative with respect to the first and 
second vector variable, respectively (thus \l/2 are dxd matrices). We also 
denote 

^ = '^{x,x), ^i = '^i{x,x), ^2 = '^2{x,x). (4.27) 



Proposition 4.3. The numerical scheme (^T^^, where \1/ satisfies ^^^2^, is 
locally exact for 

Proof: The exact discretization of the linearization of equation x = F{x) is given 
by (j3.6p . The linearization of the scheme (|4.25p (at x„ = x) reads 



sn+l Tin 



5{^ii^ + ^2in+i)+5^ (4.29) 



where 5 = 5{x) and we use ()4.27p . Identifying (j4.29p with ()3.6p we get a system 
of two equations: 

(1 _ 5^2y^{l + = e^-^' , (4.30) 



(1 _ S^2)~^S^ = (e''"^' - Ij {F'y^F , (4.31) 
where F = F{x) and F' = F'{x). Equation (|4.30p implies 

6 (^^i + ^26^"^'^ = e^^^' - 1 (4.32) 

Eliminating from (|4.32p (by virtue of (|4.26p ). we get 

5^2 (e^"^' - l) + 6F' = e'^"^' - 1 (4.33) 
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which imphes (j4.3ip (note that F = Therefore, the focal exactness imposes 
only one condition for S, namely (j4.32p . □ 

Expressing S given by fl4.28p in a more symmetric way, we obtain another, 
equivalent, form of locally exact scheme f l4.25p : 

(4.34) 

All numerical schemes presented above can be considered as particular 
cases of f l4.25p . namely: 

• expUcit Euler scheme: '^{Xn,Xn+i) = F{Xn) , 

• implicit Euler scheme: ^(x„,:e„+i) = F{Xn+i) , 

• imphcit midpoint rule: ^(x„,x„+i) = F(a;+) , 

• trapezoidal rule: \E'(x„,x„+i) = | {F{Xn) + F(x„+i)) . 

Therefore, the proof of Proposition 14.31 is valid for all these cases. 

5 Locally exact discrete gradient schemes for 
one-dimensional Hamiltonian systems 

The discrete gradient scheme (or modified midpoint rule) is a conservative 
integrator which has been used since many years for simulating dynamics 
of systems of particles [HI |2T]. More recently discrete gradient methods 
have been developed in the context of geometric numerical integration [23], 
see [131 [So]- In particular, Quispel and his co-workers constructed numeri- 
cal integrators preserving integrals of motion of a given system of ordinary 
differential equations [24l |3ll [32]. In this section we consider Hamiltonian 
systems with one degree of freedom: 

X = Hp , p = -Hr, , (5.1) 

where H = H{x,p) is a given function (sufficiently smooth), subscripts de- 
note partial differentiation and the dot denotes the total derivative respect 
to t. The Hamiltonian H{x,p) is an integral of motion (the energy integral). 
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In order to make this paper more self-contained, in section 15.11 we present 
a locally exact symmetric discrete gradient scheme, first obtained in [TT] . 
In section 15.21 we derive new locally exact scheme modifying the coordinate 
increment discrete gradient method. 

5.1 Locally exact symmetric discrete gradient scheme 

Following |TT], we consider a class of non-standard (compare [26]) discrete 
gradient schemes for the system (15. ip : 

H{Xn+l,Pn+l) + H{Xn,Pn+l) " H{Xn+l,Pn) " H{Xn,Pn) 
Sn 2{pn+l-Pn) 

Pn+1 - Pn _ H{Xn,Pn+l) + H{Xn,Pn) - H{Xn+l,Pn+l) - H{Xn+l,Pn) 

(5.2) 

where 5„ is an arbitrary positive function of hn,Xn,Pn,Xn+i,Pn+i etc. (the 
time step is denoted by /i„). The subscript n indicates that 6n may depend 
on the step n. The discrete system ( 15. 2p is manifestly symmetric (time- 
reversible). 

Lemma 5.1. The scheme ( 15.^) exactly preserves the energy integral for any 

Sn, i.e., H{Xn+l,Pn+l) = H{Xn,Pn)- 

Proof: In order to obtain the energy conservation law it is enough to multiply the 
first equation by 2{pn+i — Pn) and the second equation by 2{xn+i — Xn), and to 
subtract resulting equations. □ 

Proposition 5.2. The discrete gradient scheme ( 15. S\} with 

2 hnUJn 

— tan 

Un 2 



Sn = — tan , ujn = J H^xHpp - H^^ , (5.3) 



(where Un is evaluated at x,p) is locally exact. 

Proof: We linearize ()5.ip . substituting x = x + S,, p = p + i]- 

i = Hp + Hp^e. + HppT] , T] = -H^ - H^^^ - H^pT] . (5.4) 
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The exact discretization of the system (j5.4p is given by 



^"+1 ) = e^-^' ( M + (e^-^' - l) {F'r'F , (5.5) 



n 

(compare Proposition [27TJ, where 



F' = ( ^""^ ) • (5-6) 



Then, we hnearize the system (|5.2p around x,p, obtaining 

— = Hp + -Hxp {in + in+l) + 2^PP + ^n+l) , 

— = —Hx — -Hxx {in + in+l) — TlHxp {Vn + Vn+l) 

On 2, 2 



(5.7) 



where x„ = x + ^n, Pn= P + Vn and partial derivatives Hx-, Hp, Hxx,Hxp and Hpp 
are evaluated at x,p. The system ()5.7p can be rewritten in the matrix form 

where F and F' are defined by (j5.6p . We easily verify that 

{F'f = -u;l, u:l = HxxHpp-H%. (5.9) 

where here (and in many other places) we omit the unit matrix (i.e., we write 
instead of w^/, etc.). Comparing (jS.Sp and (jS.Sp we obtain local exactness 
conditions: 

e'^n^' = (1 - i5„F')-Hl + I'^n.n , 

(5.10) 

(g/^nF' _ i)^F')-^F = (1 - \5nF')~HnF . 

Substituting the first equation into the second one, we get an identity. Therefore, 
it is enough to consider the first equation: 

eh^F' _ UnF'e^-^' = 1 + UnF' =^ 5n = 2(F')"' tanh ^/i„F' . 

Therefore, 5n depends analytically on {F'Y and, by virtue of (|5.9p . 5n is propor- 
tional to the unit matrix. Hence (5„ is indeed a scalar function, given by 

1 2 1 

6n = 2{F')~^ tanh -/i„F' = — tan -/i„a;„ , 

YJ 



which ends the proof. □ 

As usual, we may take either x = x^, p = (to keep the scheme 
symmetric), or x = x„, p = Pn (to minimize the computational cost). 

We point out that the formula (15. 3p implies some limitations on hn in the 
case Un € IR. Certainly we have to require hnOOn 7^ vr + 2nM (M e N), or 
even hnUJn < tt. The last inequality is quite reasonable because it means that 
hn < where T„ = 271 / Un is a. corresponding period, compare [9]. 

5.2 Locally exact coordinate increment discrete gradi- 
ent scheme 

In the previous section we used symmetric form of the discrete gradient. The 
case of coordinate increment gradient (see [2D]), although of simpler form, 
turned out to be more difficult in the context of locally exact modifications. 
However, we succeeded to derive such modification also in this case. 

We consider the following non-standard numerical integrator for the Ha- 
miltonian system (15. ip : 

H{Xn+l,Pn+l) - H{Xn+l,Pn) 



Pn+1 Pn j^j^n 
Pn+1 -Pn _ H{Xn,Pn) " H{ 

5 Pn ) 

-^n+l -^n 

where, similarly as in the formula (15. 2p . 5n is an arbitrary function. 
Lemma 5.3. The scheme Ii5.11\) exactly preserves the energy integral for any 

5n, i.e., H{Xn+l,Pn+l) = H{Xn,Pn)- 

Proof: We multiply the first equation by Pn+i — Pn and the second equation by 
Xn+i — Xn- Then, we subtract resulting equations. □ 

Proposition 5.4. The scheme ( (5. iij) is locally exact for 6n given by 



2 



~ Uncot^ + H^, ' ^rr-^H^.H,,-H^^^, (5.12) 



where derivatives of H are evaluated at x,p. 
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Proof: Using notation from section \5A\ we linearize (j5.1ip around {x,p): 

S,n+1 — £,n = ^nHp + \6nHpp{rin + ??n+l) + Hxp^n+l , 

which can be rewritten as (compare (j5.6p ): 

(5.14) 

Requiring that (|5.14|) is identical with the exact discretization (|5.5|) we get: 



(5.13) 



{l - ^Hxp6n - ^6nF'^ e'^"^' = 1 - ^Hxp6n + ^dnF' , 

1 - ^HxpSn - ^5nF'^ (e'^"^' - l) (F'r'F = 5nF . 



(5.15) 



Substituting the left-hand side of the first equation into the second equation, we 
get an identity. We use the first equation to compute (5„: 

6n = 2 (e^-^' - l) ({e^-^' + l)F' + {e^-^' - l)Hxp) ' , (5.16) 



which reduces to 



5n = 2{Hxp + F'cot^) . (5.17) 



Finally, we observe that function f{z) = zcotz is even and depends analytically 
on z^. Moreover, [F')"^ = —uj'^ is proportional to the unit matrix. Therefore 6n is 
a scalar function (which is compatible with our assumption), given by (I5.12p . □ 



5.3 One-dimensional separable Hamiltonian systems 

As a simple example we consider the case defined by H^p = 0, i.e, 

H = T{p) + V{x) . (5.18) 

In this case the coordinate increment discrete gradient becomes identical with 
the symmetric discrete gradient and the scheme (15. 2 p reduces to 

Xn+l - Xn _ T{pn+l) - T{pn) 
8n Pn+l - Pn 

(5.19) 

Pn+1 - Pn _ V{Xn+l) - VjXn) 
•^71+ 1 -^n 
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Local exactness yields 



5„ = /i„tanc^ , Un = VT"{p)V"iq) . (5.20) 

Separable Hamiltonians flS.lSp are associated with many mechanical systems 
with one degree of freedom. The case T{p) = was considered in previous 
papers [H |9l [TT] , where many numerical experiments were reported. Assum- 
ing X = x„, p = p„ we get a scheme called GR-LEX, while x = \ (x„ + x„+i), 
P = \ iPn + Vn+i) yields GR-SLEX [TT]. The system f lS.ip is symmetric (time- 
reversible). GR-SLEX preserves this property, while GR-LEX does not. 

The discrete gradient schemes GR and MOD-GR are of second order. 
Locally exact discrete gradient schemes have higher order: GR-LEX is of 
3rd order and GR-SLEX is of 4th order, see P|. Numerical experiments 
presented in [2l|TT] have shown that the accuracy of GR-LEX and GR-SLEX 
is higher by several orders of magnitude when compared with the standard 
discrete gradient method (while the computational cost is higher at most 
several times, usually much less). Locally exact modifications turn out to be 
of considerable advantage. 



6 Locally exact discrete gradient schemes for 
multidimensional Hamiltonian systems 

This section contains main results. We extend results of section construct- 
ing energy-preserving locally exact discrete gradient schemes for arbitrary 
multidimensional Hamiltonian systems in canonical coordinates: 

dH dH 

X = TTT , P = ■ 6.1 

op'^ ox" 

where k = 1, . . . ,m. We obtain two different numerical schemes using either 
symmetric discrete gradient or coordinate increment discrete gradient. 
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6.1 Linearization of Hamiltonian systems 

We denote x = (x\ . . . , x'")^, p = {p^, ■ ■ ■ ,p"^), etc. The linearization of 
(16. ip around x,p is given by: 



k=l k=l 
m m 

k 

k=l k=l 



(6.2) 



or, m a matrix notation 



^ = Hp + Hpx^ + HppT) , 

(6.3) 

V = ~Hx — Hxxi — HxpT] , 

where derivatives of H are evaluated at x,p. Note that mxm matrices H^x, 
Hxpi Hpxi Hpp satisfy 

= Hpp , = Hxx , H^p = Hpx (6-4) 



Equations (16. 3 p can be rewritten also as 

dt\V J \V 
where 

\ ^2 / \ Hxp J 

Corollary 6.1. The exact discretization of linearized Hamiltonian equations 
lid. 5\) is given by 



^n+l \ _ ^h„F' f \ . ( ^hnF' i \ 



J \ Vn 



Proof: The exact discretization of (|6.5p is given by ()6.7p which follows immediately 
from Corollary 13.31 □ 
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6.2 Discrete gradients in the multidimensional case 

Considering multidimensional Hamiltonian systems we will denote 

where x E M™, p E M™, y E M^™, etc. In other words, 



A discrete gradient, denoted by V-ff (?/„,?/„ 



+1^ 



V A?' A^' ■ ■ ■ ' J = V A^' A^ ) 



is defined as an R -valued function of Vn^Vn+i such that 

t^i (6.10) 
VHiy,y) = iH,,Hp) , 
where i^i. Hp are evaluated at y = (:r,p) and we define 

WH{y,y) = \imVH{y,y) (6.11) 

when necessary. 

Discrete gradients are non- unique, compare [131 EHl [21] ■ Here we confine 
ourselves to the simplest form of the discrete gradient, namely coordinate 
increment discrete gradient [20] and to its symmetrization. The coordinate 
increment discrete gradient is defined by: 



AH H{yl^^, yl yl . . . , yl"") - H{yl yl yl...,y 



2m\ 
5 y-n J 



Ay^ yl+i - yi 

AH _ H{yl^„ yl^„ yj . . . , y^) - Hjy'^^,, yl . . . , y^) 

Vn+l Vn -|^2) 



AH _H{yl_^„y:^ 



n+l' 



Ay 



2m 



n,2m „,2m 
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where, to fix our attention, y is defined by fl6.8p . In fact, we may identify with 
y any permutation of 2m components x^,p^ . Thus we have (2m)! discrete 
gradients of this type (in particular cases some of them may be identical). 

Having any discrete gradient we can easily obtain the related symmetric 
discrete gradient 

VsH{yn,yn+i) = \ {VH{yr„yn+i) + VHiyn+i^y^a)) ■ (6.13) 

One can easily verify that VsH satisfies conditions (16.101) provided that they 
are satisfied by VH. 



6.3 Linearization of discrete gradients 



In order to construct locally exact modifications we need to linearize discrete 
gradients defined by (I6.12p and f l6.13p . 

Lemma 6.2. Linearization of the coordinate increment discrete gradient 
^UTTE) around y yields 



VH{yn,yn 



1 



■ 2 

where we denoted Vn = Vn — §, and 



(6.14) 



/ iff 

TT 



A 





1 TT 

2 y'^y^ 













\ 



H y2m—l yl 



Hy2m-ly2 
Hy2my2 



y2m — l y2m — l 
H y2m y2m—l 



(6.15) 



B 







y^y 

-H 2 2 

2 y y 












H yly2m—l 

Hy2y2m-1 



Hyly2m 

Hy2y2m 



"2 Hy2m ~ 1 y2m — 1 



H y2in-~l y2m 
'^Iiy2my2m 



Proof: We denote yi = (y^+i, • • • , y^+i, yn'^\ . . . , y'^Y ■ In particular, y° = y„ 
and = y„,+i. Expanding H{yii) around yi^^ , we get 

Hiyi) = F(yr')+i^j;.(yr')(y^+i-y^;)+^^y.j,.(yr')(y^+i-yD'+... (e.ie) 
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Hence 

= Hiyi)-H{yi-') ^ (.^-1) ^ I iyi''){yU, - y^) + . . . (6.17) 
Then, from the definition of we have 



(6.18) 



and, taking it into account, we rewrite (|6.17|) as 

i-i 2m ^ 



k=l k=j 



(6.19) 

which is equivalent to (|6.14p . (|6.15p . □ 
Lemma 6.3. Linearization of the symmetric discrete gradient yields 

VsH{yn,Vn+l) ^ Hy + ^Hyy {Un + Un+l) , (6.20) 

where Hyy is the Hessian matrix of H , evaluated at y. 

Proof: We observe that A + B = Hyy and then we use (j6.13p and (j6.14p . □ 

6.4 Conservative properties of modified discrete gra- 
dients 

In the one-dimensional case the corresponding locally exact modification is 
clearly energy-preserving, compare section 15.11 In the general case, conser- 
vative properties are less obvious. In this section we present several useful 
results. 

Lemma 6.4. We assume that a 2m x 2m matrix A (depending on h and, 
possibly, on other variables) is skew-symmetric (i.e., A'^ = —A) and 

lim ^ = S , 5 = 1 T n ) • (6-21) 
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Then, the numerical scheme 

yn+i-yn = A^H , (6.22) 

where Pn is defined by Ii6. 8\) and VH satisfies Ii6.10\) . is a consistent integrator 
for 1^6. 1\) preserving the energy integral up to round-off error. 

Proof: The consistency follows immediately form (|6.2ip . The energy preservation 
can be shown in the standard way. Using the standard scalar product in M?"^, we 
multiply both sides of HiG^ by VH 

{VH I y„+i - y„) = {VH \ AVH) . (6.23) 

By virtue of (|6.1U|) the left-hand side equals H{yn+i) — H{yn). The right hand 
side vanishes due to the skew symmetry of A. Hence H{yn+i) = H{yn). □ 

Lemma 6.5. We assume that 2m x 2m matrix 6 is of the following form 
S —a \ T T ^ 



*=(p F } • " " = l C-^^' 

(where 5, a , p are m x m matrices) and VH is (any) discrete gradient. 
Then, the numerical scheme given by 

Vn+i -yn = eSVH (6.25) 

preserves the energy integral exactly, i.e., H{xn+i,yn+i) = H{xn,yn)- 

Proof: We observe that dS is skew-symmetric, and then we use Lemma |6.4[ □ 

We easily see that fl6.25p reduces to the standard discrete gradient scheme 
when we take 9 = (i.e., 9 is proportional to the unit matrix). 



Lemma 6.6. // 6 is of the form i6.24\) and z i— )■ f{z) is any analytic 



function, then f[6) is also of the form i6.24\} 



Proof: First, we will show that the conditions (j6.24p are equivalent to 

e'^ = s-^es . (6.26) 
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Indeed, assuming a general form of 6, e.g., ^ ~ ^ ^ 1 ^ that the con- 

straint (j6.26p is equivalent to 7 = 5^, p"^ = — p and a'^ = —o. Then the proof is 

00 

straightforward. Assuming j{z) = ^^a^z", we obtain 

fe=i 

(CXD \ ^ 00 / 00 \ 

ake' = ^ a, (S-'eS) ' = 5-M ^ aj' S , (6.27) 
k=l / k=l \k=l / 

i.e., fief = s-^fie)s. □ 



Corollary 6.7. // 6'^ = S ^6S and f is an analytic function, then the 
scheme t/n+i —yn = f{'9)S'VH preserves exactly the energy integral H. 

6.5 Locally exact symmetric discrete gradient scheme 

We begin with the symmetric case because section |5] suggests that the coordi- 
nate increment discrete gradient case is more difficult. In the symmetric case 
a locally exact modification is derived similarly as in the one- dimensional 
case. 

Proposition 6.8. The following modification of the symmetric discrete gra- 
dient scheme is locally exact at y: 

yn+i-yn = 0nSWsH , (6.28) 

where 

e„ = 2(F')-4anh^ , (6.29) 
and F' , given by ( fg. 6\) . is evaluated at y = y. 

Proof: We are going to derive (j6.29p . assuming that 9^ depends on y and h. By 
virtue of Lemma 16.31 the linearization of ()6.28p is given by 



^n+l — — OnS ^Hy + -Hyy{Vn+\ + t'n)^ • 



(6.30) 
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Taking into account (j6.6p we transform (j6.30p into 

(^1 - l^nF'^ =[l + \6nF'^ Vn + QnF . (6.31) 

The scheme (j6.28p is locally exact iff (j6.3ip coincides with (j6.7|) . Therefore, we 
require that 



(6.32) 



1 - i^„F') (e^-^' - l) {F'r^F = dnF . (6.33) 

From equation (j6.32p we can compute On which yields (j6.29p . Equation ()6.33p is 
automatically satisfied provided that (|6.32p holds. □ 

Proposition 6.9. The numerical scheme Ii6.28\) with On given by ( fg. 29^) is 
energy-preserving. 

Proof: We have F' = SHyy. Therefore, {F')'^ = -HyyS = -S-^F'S, and, as a 
consequence 

{{F')Y = S-^F'fS . (6.34) 

It means that {F')'^ has the form (j6.24p . The formula (16.29P expresses 6n as an 
analytic function of {F')"^. Finally, we use Lemma |6.6[ □ 

In the one- dimensional case {F'Y is proportional to the unit matrix which 
essentially simplifies arguments presented in this section. 

6.6 Locally exact coordinate increment discrete gradi- 
ent scheme 

The symmetric form of the discrete gradient leads to a simple form of the 
locally exact modification. It turns out, however, that starting from the 
simplest form of the discrete gradient, namely coordinate increment discrete 
gradient pO] , we also are able to derive the corresponding locally exact mod- 
ification. 
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Proposition 6.10. The following modification of the coordinate increment 
discrete gradient scheme is locally exact at y: 



Vn+l -Vn = OnSVH , 

where VH is given by 1^6.12) . 



9n = 2\ SR + F'coth 



F' is given by h6.(^) (i.e., F' = SHyy), and, finally R = A — B , i.e., 

( Hyly2 . . . Hyly2m-1 Hyly2m ^ 

Hy2yl ... Hy2y2m~l H y2 y2m. 

R = 



(6.35) 



(6.36) 






-H 1 2 


Hyly2m. — 1 


Hyly2in 


H 2 1 





— Hy2y2m~l 


— Hy2y2m 


Hy2m~lyl 


Hy2m-ly2 





Hy2m — 1 y2m. 


H y2m yl 


Hy2my'2 


Hy2my2m—\ 






(6.37) 



F' and R are evaluated at y = y. 



Proof: We are going to derive (j6.36p . assuming that On depends on y and h. By 
virtue of Lemma 16.21 the linearization of (j6.35p is given by 



l/„+l -Un= dnS{AUn+l + BUn) + dnSHy 

Hence, taking into account that SHy = F, 

(1 - dnSA) Vn+l = (1 + BnSB) I/„ + dnF . 



(6.38) 



(6.39) 



The scheme ()6.35p is locally exact iff (|6.39p coincides with (|6.7p . Therefore, we 
require that 



(1 - dnSA) e 



h„F' 



1 + enSB , 



(1 - enSA) (e'^"^' - ij {F'r'F = OnF . 
Inserting (j6.40p into ()6.4ip we get 
enS{B + A){F')-^F = enF , 



(6.40) 
(6.41) 

(6.42) 
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which is identicahy satisfied by virtue oi A + B = Hyy = S ^F' , compare (j6.15p . 
The remaining equation, ()6.40p . defines On'. 

On {sAe^-^' + SB^ = e^"^' - 1 . (6.43) 

In order to simplify (j6.43p we introduce R = A — B {R is antisymmetric because 
B = A^). Then, taking into account SA + SB = F', we get 

SA = -F' + -SR , SB = -F' - -SR . (6.44) 

2 2' 22 ^ ^ 

Substituting it into (j6.43p we complete the proof. □ 

Proposition 6.11. The numerical scheme Ii6. 35\) with On given by lid. 3b]) is 
energy-preserving, i.e., H{yn+i) = iJ(|/„). 

Proof: We have 

C = 2 (^{SRf + (^F' coth = S-^dnS , (6.45) 
because 

{SRf = {-R){-S) = S-^ (SR) S (6.46) 
and, by virtue of Lemma 16.61 



F' coth ^) ^ = I^F' coth ^) S , (6.47) 
where we took into account (j6.34p . Then, we use Corollarv 16.71 □ 



6.7 Separable Hamiltonians. Multidimensional case 

In the case of one degree of freedom the assumption Hxp = simplifies final 
formulas and yields the same results for both considered gradient schemes, 
see section [5731 In the multidimensional separable case, i.e., 

H{x,p) =T{p) + V{x) , (6.48) 

a considerable simplification occurs only for the symmetric discrete gradient. 
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Proposition 6.12. The numerical scheme 

-X„) = V sT{Pn,Pn+l) 

(6.49) 

preserves exactly the energy integral (for anymxm matrix 6n), i-e., T(j9„) + 
V{Xn) does not depend on n. This scheme is locally exact for 

Sn = 2n-' tan , nl = Tpp{x,p)V,,ix,p) . (6.50) 



Proof: The first part of the Proposition follows directly from Lemma 16.51 (in this 
case a = p = 0, and 6 = Sn)- The second part is a consequence of Proposition 16.81 
We have F = {Tp, -V^Y and 

where all quantities are evaluated at {x,p). We denote = Tj^Vxx, and then 

iF'f ^ - ( f ) (6.52) 

and 

. ( ( ^^^^^ ) ^ (..3, 

We complete the proof applying Proposition 16.81 compare also Proposition 16.91 □ 



Formulas of Corollary 12.41 are strikingly similar to those given in Proposi- 
tion [61121 This is due to the fact that locally exact discretizations applied to 
linear systems yield exact integrators. Indeed, the harmonic oscillator (12.61) 
is a special case of (I6.48P for 

Tip) = ^ = l(a; I n'^x) - (x I o) (6.54) 

Then Tpp = 1 and Vxx = f^^- Hence fl'^ = Q and S'^ = 5n- Proposition 12.51 is 
an obvious consequence of Proposition 16.121 
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7 Concluding remarks 



We presented a new construction of very accurate numerical schemes based 
on the notion of local exactness. This notion is known (although under 
different names) since almost fifty years. The original application, see |29] . 
has been confined to the exact discretization of linearized equations, compare 
section 13.31 We obtain in this way one particular locally exact integrator, 
which has some advantages (e.g., it can serve as a good predictor, (|9], section 
V.A) but lacks geometric properties and stability of, for instance, locally 
exact discrete gradient methods [9]. 

Our approach has two new features. First, we modify known numerical 
schemes in a locally exact way. Any numerical scheme admits at least one 
(usually more) natural locally exact modification, see section HI Second, we 
try to preserve geometric properties of the original numerical scheme. This 
task is not trivial. In this paper we present one successful application: lo- 
cally exact modifications of discrete gradient methods for canonical Hamilton 
equations. We are able to preserve exactly the energy integral and, in the 
same time, increase the accuracy by many orders of magnitude. Another 
advantage is a variable time step. Unlike symplectic methods (which work 
mostly for the constant time step) discrete gradient methods admit con- 
servative modifications with variable time step. Therefore, one may easily 
implement any variable step method in order to obtain further improvement. 

We devoted a lot of attention to the one- dimensional case, presenting it 
independently from the general, multidimensional case. The reason is not 
only pedagogical but also practical. One-dimensional case already proved to 
be successful [H |TT]. The proposed modification, although more expensive 
(but only by a dozen or so percents), turns out to be more accurate even 
by 8 orders of magnitude in comparison to the standard discrete gradient 
scheme. Therefore, in the case of one degree of freedom our modifications are 
very efficient. In multidimensional cases the relative cost of our algorithm is 
higher, but still we hope that our method will be of advantage. Modifications 
proposed in this paper contain exponentials of variable matrices. Similar 
time-consuming evaluations are characteristic for all exponential integrators 
and in this context effective methods of computing matrix exponentials have 
been recently developed [T6t [28]. 

We presented locally exact modifications from a unified theoretical per- 
spective. There are many possible further developments. First of all, we plan 
to apply our approach to chosen multidimensional problems, testing the ac- 



31 



curacy of locally exact schemes by numerical experiments. Then, we would 
like to extend the range of applications. Throughout this paper we assumed 
the autonomous case, x = F{x). The extension on the non-autonomous case 
can be done along lines indicated already in the Pope's paper |29]. A separate 
problem is to obtain in this case any locally exact modification with geometric 
properties. Another open problem is the construction of locally exact (or, at 
least, linearization-preserving) deformations of generalized discrete gradient 
algorithms preserving all first integrals (see [21] )• Linearization-preserving 
integrators form an important subclass of locally exact schemes. It would be 
worthwhile to study linearization-preserving modifications of geometric nu- 
merical integrators, especially in those cases when locally exact are difficult 
or impossible to construct. 
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